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ABSTRACT 


The scope of the VIP3D viscous/potential flow analysis 
method has been extended to include part-span, high-lift devices 
on general planforms. The description Df the modified code in- 
cludes details of a doublet subpanel technique in which panels 
that are close to a velocity calculation point are replaced by 
a subpanel set. This treatment gives the effect of a higher 
panel density without increasing the number of unknowns. In 
particular, the technique removes the close-approach problem of 
the earlier singularity model in which distortions occur in the 
detailed pressure calculation near panel corners . Removal of 
this problem allowed a complete wake relaxation and roll-up 
iterative procedure to be installed in the code. The geometry 
package developed for the new technique and also for the more 
general configurations is based on a multiple patch scheme. 

Each patch "las a regular array of panels, but arbitrary relation- 
ships are allowed between neighboring panels at the edges of 
adjacent patches — this provides great versatility for treating 
general configurations. Preliminary tests of the modified code 
are encouraging, but further tests are needed, particularly for 
configurations with part-span, high-lift devices. 
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SUMMARY 


The scope of the VIP3D program for the viscous/potential 
flow analysis of multi-element wings has been extended to cover 
more representative high-lift configurations. Particular ob- 
jectives included part-span, high-lift devices and an improved 
representation of the interference between wakes and surfaces. 

To fulfill these objectives, each wake has to be fully relaxed 
(i.e., made force free) and allowed to roll up; however, this 
gives rise to close-approach problems associated with the poten- 
tial flow panel methods. The close-approach problem arises be- 
cause detailed pressure calculations near a surface represented 
by singularity panels have distortions near the panel edges. 

Since these distortions have a serious impact on free-wake analy- 
sis, ways of removing the close approach problem were investi- 
gated in the two-dimensional flow case. The investigation lead 
to the development of a doublet subpanel technique in which 
panels close to the velocity calculation point are represented 
by a set of subpanels. Subpanels are generated on the inter- 
polated surface and have singularity values interpolated from 
local panel values. Subpanels, therefore, have the effect of 
higher panel density without increasing the number of unknowns. 

The three-dimensional form of the doublet subpanel tech- 
nique and also the objectives for part-span, high-lift devices 
on general planforms required the development of several new 
geometry routines. A versatile and user-oriented geometry 
package was developed based on multiple patches of panels. Each 
patch has a regular array of panels but arbitrary relationships 
are allowed between neighboring panels at the edges of adjacent 
patches. The patch, panel and subpanel arrangement is continued 
downstream on wakes which are allowed to relax and roll up in an 
iterative procedure in the potential flow code. 


Preliminary test cases of the modified code are encouraging 
and compare closely with earlier solutions; however, furt^ r 
tests need to be carried out, particularly for the case of part- 
span, high-lift devices. 
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1 . 0 INTRODUCTION 


At present it is possible to compute the flow over simple 
wings and wings with full~span high-lift devices up to incipient 
flow separation using a viscoua/potential flow iteration method 
(Ref. 1) . The computer code, VXP3D, which performs this analy- 
sis, was developed by Analytical Methods, Inc. while under con- 
tracfi' to the large-scale Aerodynamics Branch at the NASA Ames 
Researcii Center. The program is designed as a special purpose 
tool for the analysis of high-lift wing configurations. 

Under a follow-on to the contract, the potential flow part 
of the computer code has been modified to allow the method to 
be applied to problems that are more representative of the high- 
lift configurations encountered on today's aircraft. The pri- 
mary objectives were to allow ; ■ ► t-span high-lift devices and 
to improve the modeling of close interference, such as exists 
between a surface and a vortex wake. 

The modified potential flow code and preceding exploratory 
work are described in this report. The new objectives required 
the removal of the close-approach problem associated with prac- 
tical potential flow codes; this problem and ways of treating 
it are examined in Section 2.0. A doublet subpanel technique 
which resulted from the investigation in two-dimensional flow 
is described in Section 3.0, together with some test cases. The 
extension of this technique for the VIP3D code and the modified 
geometry routines are discussed in general terms in Section 4,0, 
while Section 5.0 includes details of the geometry package. The 
latter has undergone extensive modification to conform with the 
new singularity model and also to deal with the more general con- 
figurations. Section 6.0 describes the treatment of the surface 
doublet distribution. Some preliminary tests of the modified 
parts of the code are described in Section 7.0 The general lay- 
out of the VIP3D program, and in particular the viscous routines, 
remain as described in Reference 1. 


‘^'Contract NAS2-8788. 
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2.0 THE CLOSE-APPROACH PROBLEM AND ITS TREATMENT 


2 . 1 Occurrence of the Cloae-Approach P roblem 

Viscous effects and the requirement of a force-free wake 
introduce strong non-linear effects into the calculation of the 
flow about high-lift configurations. The method of Reference 1 
is a practical approach to calculating such flows using a viscous/ 
potential flow iteration, but the force- free wake requirements are 
only partially satisfied since the wake relaxation is limited to 
vertical displacements. 

Tne new objectives for more general configurations and 
part-span high-lift devices make it imperative that the full 
wake relaxation and roll-up be included (Refs. 2 and ") • It is 
important not only to obtain the correct waltn Location relative 
to downstream components in a mu.lti-eloment configuration, but 
also to get correct force-free orientation as each wake leaves 
a wake-shedding element carrying moderate to high lift levels. 

If the wake is not force-free as it leaves the surface, then 
surface pressures can be affected; for example, a cross-over 
between upper and lower surface pressures can bo calculated 
near the wing tip trailing edge when using a flat, chordwise 
wake. Figure 1 (from Ref. 4) shows an extreme example of this 
crossover effect with a fixed wake in a high-lift calculation 
and emphasizes the need for a fully relaxed wake. 

Fully relaxed wake calculations on high-lift configurations 
can lead to situations where free vortex sheets pass very close 
to the lifting surfaces, and thi.s can result in a breakdown of 
the surface singularity method. This CLOSE-APPROACH problem 
was the reason behind the restriction to vertical relaxation in 
the method of Reference 1. 

The failure of surface singularity methods in a close- 
approach situation arises because of practical reasons: the 

airfoil geometry and flow distribution have to be approximated 
and the boundary conditions can bo applied only at a finite 
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number of points. Between the boundary condit’on points, the 
flow is unconstrained, and so the detailed flow distribution in 

the model generally has local distortions. These flow distor- 
tions can influence the relaxed-wake calculations i’n two ways; 

(i) Flow distortions associated with the free vortex sheets 

« 

can influence the boundary conditions on the surface, and 
hence affect the singularity solution. 

(ii) Flow distortions associated with the fixed surfaces can 
influence the calculated location of the relaxed wake — 
sometimes even causing penetration at the fixed surfaces. 

The severity of the problem increases with circulation 
level, and is particularly bad in regions of high curvature and 
high pressure gradient. The close-approach problem — which is 
present to a varying degree in all existing surface singularity 
methods — can be alleviated (usually at the expense of more com- 
puting effort) by increasing the number of control points (i.e., 
panels) and/or going to higher-order models. 

The extent of the close-approach region (i.e., the *gion 
in which flow distortions are significant) was investigated in 
Reference 5 for discrete singularities. It was established that 
the close-approach prob.lem persists to a distance' of about one 
panel size away from the discretized sheet. An extension of 
that investigation was conducted in the present work for the 
case of a distributed singularity model; namely, the symmetrical 
singularity method (Ref. 6) , employing piecewise linear vorti- 
city and a constant source on flat panels. The symmetrical 
singularity concept minimizes singularity strengths (Ref. 6) , 
and hence minimizes the flow distortions; even so, large pres- 
sure deviations wore calculated in a detailed analytic survey 
of surface pressures near the leading edge of a GA(W)-1 airfoil 
at 20*'’' incidence. Figure 2(a). The location of panel edges is 
identified and a smooth line passing through control point 
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values emphasizes the extent of the region in which the close- 
approach problem is significant for this model. 

Figure 2 (b) shows two additional pressure surveys taken 
above and parallel to the surface at heights corresponding to 
approximately a quarter- and a half-panel size. As the pressure 
survey moves away from the surface, the extreme distortions at 
the panel edges quickly reduce to small "bumps" in the detailed 
distribution. (With increasing distance from the inclined sur- 
face, these bumps appear to . e sideways in the x-wise plot 
because the corner disturbs c propagates normal to the surface.) 
Bearing in mind that the situation shown is particularly ex- 
treme (the real flow would have separated), it would seem feasible 
for the close-approach problem area to be held within a half 
panel size for this particular model. Even so, some local treat- 
ment would be necessary to restore a smooth pressure distribu- 
tion close to the surface. 

An investigation of ways of removing the close-approach 
problem is outlined in the following svibsections . This investi- 
gation considered alternative singularity models in addition to 
various close-approach techniques applicable to existing models, 
because the symmetrical singularity panel model considered above 
is not directly extendable to the three-dimensional case without 
some deterioration in its close-approach characteristics. Al- 
though most of the evaluations were performed in the two-dimen- 
sional case, the form of each model for three-dimensional flow 
was the major consideration throughout the investigation. 

2. 2 Singularity Model 

Good close-approach properties in a surface singularity 
model require continuous and smooth representations of the sur- 
face geometry and flow distribution. These requirements con- 
flict with the properties of practical singularity models, most 
of which are based on panels over which a piecewise analytic 
integration is performed to evaluate influen-s coefficients. 
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Panels are usually flat, but higher-order effects due to surface 
curvature and smooth singularity distributions may be approxi- 
mated using a series expansion in the integrand. Generally the 
close-approach problem is present to a varying degree in all 
these methods making it necessary to couple them with some 
close-approach technique to accomplish the present objectives. 

An alternative to piecewise analytic integration uses a 
numerical integration scheme based on a continuous singularity 
distribution over the curved surface. Such a method has no 
close-approach problem concerning the representation of the 
surface geometry and singularity distributions, but it does have 
a problem in evaluating the local singularity contribution when 
calculating velocities close to the surface. 

To examine a numerical in+-.' iration approach a pilot code 
was assembled based on a surface doublet distribution. A 
smooth doublet distribution is the most convenient singularity 
model for the three-dimensional lifting case because it auto- 
matically satisfies the requiremnt of continuity in the vorti- 
city distribution (vorticity being the gradient of the doublet 
distribution) . The doublet distribution and the surface geo- 
metry were represented using the biquadratic expression (Appen- 
dix A) . The numerical integration scheme used the Romberg 
method coupled with Richardson's extrapolation technique (Appen- 
dix B) . It will bo recalled that this method is based on the 
trapezoidal rule integration, but, by continuously doubling the 
number of strips, the error value can be controlled and the com- 
puting effort minimized. 

The main problem associated with numerical integration over 
a doublet sheet is the evaluat-> n of the local contribution 
near the velocity calculation point. It arises from a unique 
behavior of the doublet integrand which first goes to a large 
positive value and then to a large negative value as the velo- 
city calculation point is approached, see Appendix B. To make 
the numerical scheme more effective in the pilot code, various 
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transformations were applied in a number of regions based on 
the behavior of the integrand (Appendix B) . The Romberg inte~ 
gration scheme became very effective with this treatment and 
gave accurate solutions of the two-dimensional airfoil problem 
with execution times comparable with piecewise analytic methods. 
A three-dimensional version, however, was not as effective. Al- 
though exact agreement was obtained in comparison with a method 
using linear vorticity panels for several velocity scan lines 
approaching the surface, the overall computing time for the 
numerical scheme was an order of magnitude higher than for the 
analytic scheme. Even though it seemed feasible to reduce com- 
puting effort further by "fine tuning" the transformation 
schemes, it was clear that considerable development time would 
be needed to make the numerica] r'-heme a prortncal and fail- 
safe method for the general three-dimensional configurations 
envisaged. Effort was thereafter concentrated on investigating 
close-approach techniques that are applicable to existing singu- 
larity models. The various techniques considered arc described 
in the next subsection. 

2 . 3 C lose-Approach Techniques 

There are several close-approach techniques which can bt 
applied to existing singularity models to obtain smooth pres- 
sure calculations between control points. Techniques considered 
in the present work are described below. Details of some of 
the evaluations are included for the calculated surface pres- 
sures near the leading edge of a Joukowski airfoil at 10° in- 
cidence; the overall pressure distribution calculated by the 
symmetrical singularities method (Ref. 6) is given in Figure 3. 
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Figure . Comparison oi I'n'usure Distributions on a 
Lloukowski Aiiioil vit 10'' Incidence. 


^ ^ Interpolation 

A technique that has been used in the past interpo- 
lates between the control point values to obtain intermediate 
surface pressures. Likewise, off-body pressures may be obtained 
by interpolating between "good" values calculated outside the 
close-approach region — such a technique has been used success- 
fully even with discrete singularity models. 

Clearly, interpolation is a reasonable technique to apply 
to attached flow single-airfoil problems (e.g., Figure 3), es- 
pecially if the interpolation is performed in terms of surface 
distance rather than x or z; however, in the case of multiple- 
element high-lift calculations, situations exist where inter- 
polation is not applicable. For example, there are no "good" 
values where the close-approach regions of a free wake and a 
fixed surface overlap. Higher panel density may alleviate' the 
problem, but would increase the number of unknowns and — if com- 
puter storage would allow such an increase — would be more ex- 
pensive to run. 

2.3.2 Treatment of the Panel-Edge Singularity 

The main culprit behind the local flow distortions 
is a logarithmic term which can become singular at each panel 
edge. Under certain conditions the logarithmic singularity 
terms from two adjacent panel edges cancel; e.g., when the sur- 
face slope and singularity value are continuous from panel to 
panel. For the more general case, an interpolation formula 
was examined in which the induced velocity influence coeffi- 
cients from two neighboring panels were combined linearly ac- 
cording to the location of the calculation point within the 
close- approach region. Results from calculations based on a 
limiting close-approach distance of a quarter panel size from 
the corner point are presented in Figure 4 . Although consider- 
able improvement is indicated over the basic solution, the 
results are not acceptable for the present objectives. Further 
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development of this technique is possible i but any further com- 
plication from the tiodel evaluated could get very cumbersome 
in the general three-dimensional case. 

2.3.3 Internal Singularities 

Earlier investigations (Ref. 5) of discrete singu- 
larities in two-dimensional flow lead to a "submerged singularity" 
technique (Ref. 7) in which the close-approach problem region 
was enclosed within the airfoil contour by placing the singu- 
larity sheet the appropriate distance below the surface. This 
submerged distance was minimized by applying the subvortex tech- 
nique (Ref. 5) as well. (In fact, the internal singularities 
technique is best applied in combination with another close- 
approach technique — in this way tho interior singularity surface 
can be more closely related to the airfoil surface and this aids 
accuracy.) This internal singularity technique gj.ves a very 
smooth pressure distribution (Ref. 7) and is applicable to simple 
wings in the three-dimensional case, but it may be difficult to 
apply to general three-dimensional configurations. 

2.3.4 Corner Panels 

In this technique the two panels adjacent to a close- 
approach velocity calculation are temporarily replaced by three 
panels. The middle panel of the temporary set straddles the 
corner between the two control points (Figure 5(a)), and the 
other two are the remaining halves of the replaced panels. The 
singularity values for the temporary panels are obtained by 
interpolation through the basic panel values. 

Again, this technique removed the "spikes" in the detailed 
pressure distribution. Figure 6, but the resulting distribution 
is not sufficiently smooth; kinks occur at the points where the 
corner panel scheme takes over from the basic panels. (The 
limiting close-approach distance was a quarter panel size in 
this case.) 
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Figure 6. Detailed Pressure Cilculations Near Leaing Edge of a 
.Toukowski Airoil. Symmetrical Singularities Method 
with Corner Panel Model (2.3.4). 
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It will be observed in Figure 6 (and in Figure 4 I that 
the control station values are not in close agreement with the 
exact solution. Accuracy was later restored by placing the 
control points on the airfoil surface rather than on the panel 
centers, Figure 7. The corresponding corner panel model, Figure 
5(b) was then re-a'i-aluaied, but, the detailed pre*- sure distribu- 
tion still had kinks at the changeover stations. Clearly, to 
remove these kinks would require the changeover point from one 
panel scheme to the other to be moved further away from the ve- 
locity calculation point, i.e., the limiting close-approach 
distance must be increased. This consideration lead to the more 
general technique described below and based on a number of sub- 
panels. 


2.3.5 Subpanel Technique 

In this technique, each panel which contains the 
velocity calculation point within a certain NEAR-FIELD RADIUS 
from its control point is divided into a number of subpanels. 
Figure 8. Subpanel corner points are obtained by interpolation 
on the airfoil surface and singularity values by interpolation 
through the panel values . When accumulating panel velocity con- 
tributions at a given point, immediately a panel detects the 
point is within its near-field radius a perpendicular is dropped 
from the point to the surface. The projected point becomes the 
center of a subpanel, the size of which is related to the height 
of the velocity point above the surface, Figure 8. Other sub- 
p>anels are constructed with size increasing with distance from 
the point until all the near-field panels are modified. (Note; 
this way of distributing the subpanels was later changed for 
the three-dimensional case, sec 4.2.) 

The region covered by the near-field radius is necessarily 
larger than the close-approach region, Figure 9, because the 
velocity distribution induced by the temporary subpanel set 
does not match that of the basic panel until some distance away 
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(e.g., see Figure 5 in Ref. 5). Thus, when one panel is rep- 
resented by subpanels, several neighboring panels must be 
similarly treated, otherwise significi.nt "jumps" occur in the 
calculated velocity distribution as we pass over the near-field 
boundary of each panel (e.g., as observed in Figures 4 and 6). 
These jumps can be made as small as we please by increasing the 
near-field radius. 

Using the near-field radius as the criterion for generating 
a subpanel set rather than the close-approach distance can re- 
sult in unnecessary use of subpanels. For example, the calcula- 
tion point in Figure 9 is inside the near-field radius of 
several panels yet is outside the close-approach regions. 

Clearly, the subpanel model is superfluous in such a situation, 
but individual panels would noi- ho. "certain" of that fact unless 
the geometric relationship between the point and all the pane.ls 
was tested at the beginning. Since such a test would duplica^ 
the evaluation of some of the geometric quantities in the panel 
influence coefficients calculation, it was decided to adopt the 
panel-by-panel, near-field radius test and accept the occasional 
"overkill" situation illustrated in Figure 9. 

For convenience, the near- field radius is expressed as a 
factor applied to each panel length. In this way the factor 
has one input value for the entire calculation, yet the local 
near-field radius varies according to the size of each pane}.. 

While the subpanel technique is an extension of the corner 
panel technique, it is also a higher-order form of the subvor- 
tex technique described in Reference 5. With this technique, 
velocity calculation at any arbitrary point on the surface or 
in the flow field never experiences a panel corner problem since 
the generated subpanels are always in the ideal relationship with 
the calculation point. The number of subpanels required for the 
near-field calculations appears to be very small. For example, 
in the case of a point on the surface, the total number of 
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influence coefficients evaluated (panels plus subpanels) is 
usually of the order of five more than the number of basic panels 
representing the airfoil. The subpanel pilot code was initially 
based on the symmetrical singularities model (Ref. 6); i.e., 
panels of linear vorticity and constant source. This was later 
converted to a doublet model which is described in the next 
section. 
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3.0 DOUBLET SUBPANEL TECHNIQUE 

3.1 Singularity Model 

Before evaluating the subpanel technique in morr. detail in 
the two-dimensional flow case, the singularity model was changed 
to a surface doublet distribution to facilitate extension to the 
three-dimemsional case. For the purpose of evaluating doublet 
values at subpanel centers, the doublet distribution is described 
by a biquadratic interpolation curve (Appendix A) passing through 
known values at panel control points. Figure 10(a). Thus, for 
each subpanel, the central doublet value is obtained in terms of 
four local panel doublet values, 

4 

U = X) ' X 

i=l 

The biquadratic multipliers, are evaluated at the subpanel 
center and are functions only of surface distances (Appendix A) . 

The linear vorticity influence coefficient routine for 
panels or subpanels is retained from the earlier code (see Ap- 
pendix C) , but the vorticity is now evaluated as the gradient 
of the doublet distribution with respect to surface distance. 
Clearly, the vorticity value could be obtained by differentiat- 
ing the biquadratic expression directly (Appendix A) , but such 
treatment of high-order curves is not always reliable. An al- 
ternative approach adopted here and based on a local quadratic 
curve passing through three subpanel doublet values, Figure 10 
(b) , makes use of the fact that the panel doublet values have 
been augmented by the biquadratic interpolation for subpanel 
values; that is, the vorticity value at the middle subpanel is 

Yj^j =-au/9s = (u^ ~ 

+ Mj^(S^ - S^)/Sj^/S^ (2) 
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(a) INTERPOLATION THROUGH PANEL DOUBLET VALUES TO OBTAIN SUBPANa VALUES 

Fiaure 10. Treatment of the Surface Doublet Distribution. 



LCXAL SURFACE LENGTH 

INTERVALS BETWEEN SUBPANEL CENTERS 


(b) LOCAL QUADRATIC REPRESENTATION FOR EVALUATING VORTICITY 
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Figure 10. Concluded. 


and the vorticity gradient is 


Since each subpanel doublet value (viz., u^, vij^, y^) is 
known in terms of a set of four local panel values, Eqn. (1) , 
then the linear vorticity influence coefficients (Appendix C) 
for subpanels can be immediately expressed in terms of panel 
doublet values. 

3 . 2 Boundary Condition 

With a smooth doublet distribution on (and with doublet 
axes normal to) a closed surface there are two equivalent ways 
of applying the exterior tangential flow boundary condition; 

(i) the exterior Neumann boundary condition in which the normal 
velocity component is set to zero, and (ii) the interior 
Dirichlet boundary condition in which the interior velocity po- 
tential is set to a constant (e.g., 0). 

The Neimann boundary condition when applied to a surface 
doublet or vorticity distribution is weak in the region approach- 
ing the airfoil trailing edge. This often leads to local devia- 
tions in the solution. The problem does not occur for the al- 
ternative interior boundary condition, and so from this stand- 
point, the latter boundary condition is to be preferred. In 
addition, since the Dirichlet form works with a scalar quantity 
(velocity potential) rather than a velocity vector, computing 
effort and storage requirements are minimized. However, with 
this approach, the surface velocities are obtained from the 
gradient of the velocity poteiitial and this might introduce in- 
accuracies in the general three-dimensional case. At this time, 
therefore, and because the VIP3D program is set up for it, the 
Neumann boundary condition is used here. The possibility of chang- 
ing to the interior boundary condition will be kept in mind for 
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a future modification subject to the availability of a reliable 
way of obtaining the surface velocities in the general case. 

Use of the Neumann boundary condition requires that the 
ill-conditioning be alleviated in the equations leading to the 
solution of the panel doublet values. The cause of the problem 
can be narrowed down to a certain part of the doublet distribu- 
tion; by considering a typical doublet distribution on the sur- 
face of a lifting airfoil we can identify symmetrical and anti- 
symmetrical parts with respect to surface distance from the 
trailing edge, Figure 11. The symmetrical part can be associated 
with the airfoil thickness or displacement effect, while the 
antisyrametrical part is closely related to the circulation. Be- 
cause the normal velocity component is continuous, passing 
through the doublet sheet, it ran be shown (see Apendix D) that 
the Neumann boundary condition equations for control points near 
the trailing edge are ill-conditioned only for the solution of 
the symmetrical part. It will be recalled that when the surface 
pressures from the doublet solution become inaccurate near the 
trailing edge, it is usually observed that the pressure difference 
between upper and lower surfaces (i.e., the circulation or anti- 
symmetric effect) is correct. 

In the symmetrical singularities method (Ref. 6), the anti - 
symmetrical doublet distribution is represented by a symmetrical 
vorticity distribution (being the gradient of the doublet distri- 
bution) while the symmetrical doublet distribution influence is 
represented by a symmetrical source distribution. Both the 
source and vorticity components are solved using the same number 
of boundary condition equations as there are panels by applying 
the symmetry constraints. In an earlier code (Ref. 7) , the sym- 
metrical source distribution was actually applied at the begin- 
ning as a simple function of the rate of change of thickness. 

The applied source influence when combined w’ith the onset flow 
was found to be adequate in stabilizing the doublet solution. 

In the present work, the presence of the source singularity is 
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undesirable; it is more convenient to work with the total doublet 
distribution for the purpose of shedding circulation from free 
edges, such as flap edges, wing tip, etc. To maintain a stable 
solution with the Nexamann boundary condition, therefore, an ap- 
proximate symmetrical doiablet distribution is applied (by the 
program) at the beginning based on the flow about a symmetrical 
Karman-Trefftz airfoil at zero lift. The Karman Trefftz section 
is constructed having the same trailing-edge angle and cross- 
section area as the actual airfoil. In this way the solution 
part of the doublet distribution is concerned primarily with 
the (antisymmetric) circulation component, but has small adjust- 
ments because the applied symmetrical component is not exact 
for that particular airfoil. (The solution part could be re- 
duced even further by simply rsotting the Karman-Trefftz section 
at the appropriate angle of attack — rather than zero — when eval- 
uating the applied doublet distribution.) 

3 . 3 Test Cases 

As a searching test of the subpanel technique, a vortex/ 
surface interaction calculation was chosen in which a prescribed 
vortex was positioned close to a Joukowski airfoil. The vortex 
location was x = .15c, z - .125c, and its strength was .2 tt. The 
vortex flow was combined with an a ^ 10° onset flow. Thirty 
panels were used in a cosine spacing, and the near-field radius 
factor was set to 3. 

The ability of the subpaneling scheme to provide smooth 
velocity calculations anywhere is very apparent in Figure 12(a), 
which shows calculated streamlines. The streamline calculation 
method is one developed during the course of the work (Appendix 
E) . Three starting points were selected as shown in Figure 12 
(b) . The forward point gives a streamline that on the upstream 

part passes very close to the leading edge, and in the down- 
stream part climbs over the vortex before dropping to the air- 
foil surface which it follows very closely back to the trailing 
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edge. Details of this streamline (and the second streamline) 
in the leading~edge region are given in the inset in Figure 12 
(a) . The first streamline passes very close to the surface, well 
within the spacing of the control points. The line is very 
smooth, even though the velocity calculations have been perfor- 
med at a nximber of "arbitrary" positions. The second streamline 
is clearly very close to the stagnation streamline and essenti- 
ally follows the surface with one or two minor oscillations. As 
the calculation proceeds from the starting point, this second 
streamline hits the airfoil very steeply, and yet quickly takes 
”p the surface direction, a very searching test for both the 
streamline calculation procedure and the velocity calculation 
routine. On the downstream side, this second streamline follows 
the surface back to the trailing edge. 

The third streamline forias a closed loop round the vortex 
and does several turns (total streamline length specified is 
2.5 chords) before accumulating errors eventually allow it to 
escape downstrean^ along the airfoil surface. 

The surface pressure distributior, corresponding to this 
calculation is shown in Figure 12(b). Intermediate velocity 
calculations are indicated by triangles to distinguish them 
from the basic control-point values. These additional calcula- 
tions, made possible by the subpaneling tec’'» ique, clearly de- 
fine the details of the three suction peaks and three stagnation 
points. The control point values in some of these areas would 
have been inadequate — particularly in defining the suction peak 
located beneath the vortex. Clearly, the interpolation tech- 
nique (2.3.1) would have been incapable of constructing these 
details given only the panel control point values. An exact 
solution for this type of problem has been provided by ur. V.J. 
Roscow at NASA Ames Research Center (Ref. 8) . The solution is 
obtained by a transformation technique and gave the results for 
a CLARK- Y section. 
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Figure 13(a) compares the exact and calculated streamlines 
for a vortex and sink. The sink, which is coincident with the 
vortex, provides the necessary stabilizing force to keep the 
vortex in equilibrium in the exact analysis. The sink influence 
was included in the pilot code for this calculation. The stream- 
line calculations (using the procedure described in Appendix E) 
show very close agreement with the exact lines. Figure 13(b) 
shows the close agreement between the calculated and exact pres- 
sure distributions for this case. The calculated pressures are 
at arbitrary stations (i.e., not necessarily at control points) 
and closely represent the multiple stagnation point and peak 
suction features. 

These calculations clearly demonstrate the effectiveness of 
the subpanel technique in a px’oblem situation that is pertinent 
to the multi-element high-lift configuration analysis; viz., 
vortex/surface interaction. The main features of the technique 
are summarized below. 

(i) Subpanels offer a closer representation of curved sur- 
faces and smooth singularity distributions than is 
possible with practical panel densities. 

(ii) Subpanels give the effect of higher panel density 
without increasing the number of unknowns. 

(iii) Subpanels give a "higher-order effect", yet maintain 
simple influence coefficient expressions. 

(iv) A panel's subpanel set is used only when a velocity 
calculation is performed within a small near-field 
radius from the panel's center (e.g., within three 
panel sizes av;ay) . This minimizes computing effort. 

(v) Smooth velocity calculations are obtained with reason- 
able panel density, even in the case of the vortex/ 
surface interference problem. Features of the inter- 
ference pressure field are closely represented even 
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though the control points have not been specially 
located relative to the vortex position. 


4.0 OVERVIEW OF THE VIP3D CODE MODIFICATIONS 

4.1 General 

The new work concerns only the WBOLAY part of the overlay 
structure in the VIP3D program (Ref. 1); boundary layer routines 
are unaffected. The potential flow routines — including the geo- 
metry package — have been modified and new routines added in ac- 
cordance with the new objectives and the new singularity model. 
The flow chart for the initial calculations covering the geometry 
specification and potential flow solution is shown in Figure 14. 

The program modifications and new capability are described 
in general terms in the following subsections. Details of the 
new work are given separately in later sections. 

4 . 2 Geometric Model 

The geometry routines in VIP3D have been extended to satisfy 
the new objectives concerning more general configurations and 
also to be compatible with the new close-approach singularity 
model. For these purposes — and also as a convenience to the 
user — a configuration is now broken down into more parts than 
before. In descending order of size, these parts include as- 
semblies, components, patches and wakes, panels and, finally, 
subpanels. The nature and purpose of each of these parts is 
discussed in the following paragraphs, but not necessarily in 
the same order as in the list above. 

The program may be applied to a configuration having a num- 
ber of COMPONENTS; e.g., wing, slat, flap. A component is the 
smallest part of a configuration for which integrated load and 
moment information is given in the program output. Neglecting 
at this time the effects of structural attachments such as be- 
tween a slat and wing, components are normally regarded as sepa- 
rate parts of the configuration, but this need not be so. As a 
user convenience, the surface of a configuration may be continu- 
ous from component to component, thus allowing integrated forces 
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Figure 14. Flow Chart for the Potential Flow Calculations in the 
Modified Code. o-, 







and moments to be evaluated for just part of a wing or for a 
winglet, for example. Contiguous components are collected to- 
gether into ASSEMBLIES within the program to allow the doublet 
distribution to be described continuously over the connected 
surfaces. 

As in the earlier code, the basic unit representing the 
surface of each component is the PANEL. Panels have four 
straight sides, but one of the sides may be of zero length. In 
the new code each panel may be subdivided into SUBPANELS for the 
purpose of near-field calculations. The subpanel scheme differs 
from that described earlier in 2.3.5 by having a fixed subpanel 
set formed and stored for each panel at the beginning of the 
analysis. This change was made because a complete regeneration 
of subpanels for each near-f'io.'' ^ velocity calculation, while 
practical in the two-dimensional case — the computing effort 
being relatively small — could become unacceptable in the three- 
dimensional case where the calculations are more involved. 

15 shows a section through the new subpanel scheme 
for comparison with the earlier model in Figure 8. The subpanels, 
which are generated automatically within the code, form a "square" 
array within each panel with an odd number per side. This ar- 
rangement ensures the central subpanel of a set falls in the mid- 
dle of the panel; the center point and unit normal vector of this 
subpanel are then adopted by the panel for its control point geo- 
metry. This treatment places the panel's control point very close 
to the surface rather than on the panel; the advantage of this 
was observed earlier in Figure 7 . The subpanel sides are straight 
but the corners lie on the interpolated surface represented by a 
two-way biquadratic scheme. Compared with the panel model, there- 
fore, the subpanel set represents surface curvature more closely. 
New routines installed in the code to interpolate and dif- 
ferentiate the surface doublet distribution require information 
concerning neighboring panels. Similar information is required 
when redistributing the boundary layer displacement source values 
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Figure 15. Section through Subpanel Scheme for Three-Dimen- 
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at the panel centers, and may be anticipated for the future in 
connection with surface streamline calculations. The most effi- 
cient way to locate neighboring panels is to arrange the panels 
in a "rectangular" grid of rows and columns (i.e., equal number 
of panels in each column) . Such an arrangement, however, lacks 
versatility when considering general configurations. A com- 
promise has been reached in the new code by the use of PATCHES 
of panels, each patch having a rectangular array of panels. 
Awkwardly shaped components may be represented by several patches, 
whereas simple shapes, for example, the entire main surface of 
a swept tapered wing, may require just one patch. Figure 16 
shows a typical breakdown of a high-lift configuration into com- 
ponents and patches. Parts A, B, C and D are typical patches on 
Component 2. Thie main surface*' •'f components 1 and 3 may be 
formed by single patches; additional patches mail be used to cover 
the open ends (shaded) . 

The developed (i.e., "opened out") shape of a patch should 
be roughly four-sided to keep panel shapes and distributions 
reasonably regular. This does not exclude the presence of kinks 
in any or all of the patch sides, but kink angles should not be 
large (the upper limit has not been established, but for the 
time being 60° should be regarded as a large kink angle) . One 
side or two opposite sides of a patch may be reduced to zero 
length provided the overall patch shape is reasonably regular. 
Figure 17 gives some basic guidelines for acceptable patch shapes. 

The versatility of the patch scheme is ensvired by allowing 
arbitrary relationships to occur between neighboring edge panels 
across patch joints. Figure 18. Automatic procedures (described 
later in 5.4.3) have been developed in the code to select "pre- 
ferred" neighbors from edge panels on neighboring patches. Only 
edge panels on patches within the same assembly of components 
are considered as possible neighbor candidates. (This distinction 
is necessary in the automatic procedure to avoid the possibility 


40 



(MAIN WING) 


Figui'e 16. 


Examples of Components and Patchus on a High-Lift 
Configuration . 


41 




of, say, a slat panel being selected as a preferred neighbor to 
a wing panel.) 

A number of factors contribute to the usability cf the 
patch scheme; in particular, options are provided for automatic 
paneling and also for automatic patch generation in special situ- 
ations, o.g., wing tip. Those options, which are described later 
in Section 5, are activated by setting a few integer flags in 
the input deck; they offer a very convenient mode of operation 
for a minimum of data preparation. These options are provided 
without sacrificing the basic simplicity of the input format. 

The surface doublet distribution passes continuously onto 
WAKES representing the shed vorticity sheets from each component. 
Initial wake geometries are generated within the program follow- 
ing a small amount of user input- This info unc tion is needed to 
form an initial representation of the often complicated wake ar- 
rangements seen at part-span flap cutouts and edge separations. 

Each wake is in two parts; a near- and far-wake. The near- 
wake starts at the wake shedding line on the component and ex- 
tends downstream a short distance beyond the end of the configu- 
ration. The near-wake, which is similar to a patch, has panels 
and, for close-interference calculations, subpanels (Figure 19). 

The far-wake model extends from the end of the near-wake 
back to dovrnstream infinity and is represented by semi-infinite 
vortices (i.e., a piecewise constant doublet model). Those vor- 
tices are attached to wake panel corners at the end of the near- 
wake. The far-wake requires no detailed representation since 
it is remote from the regions of interest. 

The simple wake shape iteration in the earlier code allowed 
vertical movements only. The new singularity model has lifted 
that restriction so the new wake shape iteration, which involves 
the potential flow code only, now allows full roll-up in the 
near-wake . 
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4 . 3 Singularity Model 

The piecewise constant source and vorticity model in the 
VIP3D program has been replaced by a continuous surface doublet 
distribution incorporating the three-dimensional forms of the 
subpanel close-approach technique. With the fixed subpanel sets 
of this latest technique (4.2) it is possible for a calculation 
point to approach the edge of a subpanel. The interpolation 
technique described in 2.3.1 is therefore coupled with the 
scheme using velocities calculated above local subpanel centers. 
Interpolation within the subpanel system is acceptable because 
the associated close-approach problem area is considerably 
smaller than that for the basic panel system. 

Each subpanel is associated with a two-part doublet value 
located at its center point; the applied symmetrical part is 
evaluated at the beginning (3.2) according to the local chord- 
wise geometry, while the solution part is initially described in 
terms of biquadratic interpolation multipliers applied to a local 
set of panel doublet values. These multipliers are evaluated 
from a two-way biquadratic interpolation scheme passing through 
panel doublet values (solution part) located at control points 
where the surface boundary condition is applied. The panel's 
applied doublet part is provided by the middle subpanel in its 
set. 

The total doublet value (applied plus solution) occurring 
at the trailing edge is passed onto the wake panels. The result- 
ant doublet value (i.e., between upper and lower surfaces at 
the trailing edge) is held constant in the streamwise direction 
for each column of wake panels. 

4 . 4 Solution Procedure 

The solution procedure in VIP3D remains essentially un- 
changed except that we now solve for panel doublet values. The 
matrix routine, therefore, calls a new influence coefficient 
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procedure (described later in 6.3). When the doublet solution 
is obtained, it is immediately combined with the applied doublet 
values prior to the pressure distribution calculations. 
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5.0 GEOMETRY ROUTINES 

This section describes the way the modified geometry rou- 
tines treat the various parts of the configuration. Details are 
included to emphasize the versatility of the new routines. The 
discussion is user-oriented to help in applying the new capabil- 
ity to general configurations. 

The reference coordinate system used in the original code 
is maintained here but i =! now referred to as the general coordi- 
nate system, or G.C.S., to distinguish it from two other refer- 
ence systems, namely, the component coordinate system, or C.C.S., 
and the section coordinate system, or S.C.S. These other systems 
are described in this section and are provided as a user con- 
venience for specifying the geometry of a configuration. 

5 . 1 Components 

When defining the surface geometry of a configuration, each 
component may be described in its own local coordinate system 
for convenience. This also allows components to be relocated 
at a later date with minor changes to the input deck, for ex- 
ample, a slotted flap may be moved to a different setting. 

The component specification starts with the appropriate 
transformation information which converts from the component 
coordinate system (referred to as C.C.S.) to the general coordi- 
nate system (referred to as G.C.S.). This information includes 
(i) the translation vector, (CTX, CTY, CTZ) , which is simply the 
origin of the C.C.S. expressed in the G.C.S., (ii) the scaling 
factor, and (iii) the rotation angle, 0, about a hinge line vec- 
tor, h, Figure 20. Provision is made for the user to specify 
two points on a general hinge line vector (in the C.C.S.), other- 
wise the y-axis in the C.C.S. is used. Both the scaling and the 
rotation are applied in the C.C.S. prior to the translation. 

This component transformation is performed at the end of the 
geometry input routine, i.e., after the basic geometry of the 
complete configuration is assembled. 
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5.2 Patches 


5.2.1 Convention 

We have seen (4.2) that a patch is basically a four- 
sided shape when developed or "opened out". It must always be 
regarded as such, even if one of the sides or two opposite sides 
are made zero or if some of the sides have kinks. In the follow- 
ing discussions patches will often be regarded as rectangular — 
this is purely a convenience for discussing relationships and is 
not a shape restriction. Our view of the patch will always be 
from the outside , i.e., looking onto the surface from a point in 
the flow field. 

For convenience, a patch is defined in terms of a "chordwise" 
and a "spanwise" direction, Figuio 21. These directions are 
analogous to the conventional wing layout, but, in the patch con- 
text these directions are not restricted to the x and y direc- 
tions, respectively. For example, on a patch representing the 
wing tip, the "spanwise" direction will probably be in the wing 
chordwise direction and the "chordwise" direction will be ver- 
tical. 

Patch geometry is defined using chordwise lines called 
SECTIONS. (These are described later in 5.2.2.) A set of sec- 
tions distributed spanwise across a patch defines the patch 
surface. The convention adopted here is that points defining 
a section shape proceed from top to bottom. Figure 21. (In 
the case where a patch represents the main surface of a wing, 
this convention causes the points defining each section to pro- 
ceed from the trailing edge lower surface and finish at the 
trailing edge upper surface, i.e., as in the original program.) 

In our view of the patch, the order of the sections always pro- 
ceeds in the positive spanwise direction; however, for user con- 
venience, we allow the spanwise direction to proceed either 
from left to right or vice versa, Figure 22. As we shall see 
later (5.4.3), it is important that the program distinguishes 
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between the two possibilities so we set a patch IDENT parameter 
to +1 if the spanwise direction is selected from left to right 
(i.e., root to tip in the wing convention) or -1 if the direc™ 
tion is right to left. 

For the purpose of automatically connecting panels from one 
patch to another, it is important to identify patch sides . The 
convention adopted here is that the first and last sections de~ 
fining a patch correspond to sides 1 and 3, respectively, while 
the patch top and bottom correspond to sides 4 and 2, respect- 
ively. With this convention, the order of the sides is anti - 
clockwise when IDENT is +1 and clockwise when IDENT is -1, Figure 
22. The order of the corner points follows the same sequence as 
the sides, start;.ncf with 1 at 1.^' i-.op of sic '^ T. 

Panels and subpanels take Li e same side ajid corner point 
convention as their parent patch. For convenience, the panels 
are referred to in ROWS — which run spanwise — and COLUMNS--which 
run chordwise. Figure 22. Panel arrangements within a patch are 
referred to by ROWS X COLUMNS. 

5.2.2 Sections 

Each section of a patch may be defined in its own local 
coordinate system, referred to as the section coordinate system, 
or S.C.S. The user provides the necessary information to trans- 
form from the S.C.S. into the C.C.S. at the beginning of each 
section. This transformation is performed immediately a section’s 
geometric description is complete. This transformation is sepa- 
rate from that described earlier in 5.1 in which the completa 
component is converted into the G.C.S. (at which stage the S.C.S. 
geometry is discarded) . This double transformation — both levels 
of which are optional — offers useful flexibility when preparing 
the input data. One particular advantage is that the geometric 
relationships — especially the rotations — are kept reasonably 
simple without sacrificing generality. 
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The information required to transform from the S,C,S. into 
the C.C.S, (see Figure 23) consists of; (i) the translation 
vector, (STX, STY, STZ) , which is the position vector of the 
S.C.S. origin expressed in the C.C.S. coordinates; (ii) a scaling 
factor (default InO) which is applied in the S.C.S. ; (iii) the 
rotation angle (0, degrees) about the y-axis of the S.C.S. ; and 
(iv) the angle (4>, degrees) in the C.C.S. x-y plane, between 
the projection of the S.C.S. x-axis and the C.C.S. x-axis. 

The contour line of each section is defined by a set of 
BASIC POINTS, (BX, BY, BZ) . These points may be used directly 
as panel corner points, i, e., MANUAL PANELING, in which case . 
the user must take care over the number of input points. Al- 
ternatively, an AUTOMATIC PANELING ROUTINE, or A.P.R. (described 
in 5.4.1) may be activated which interpolates through the basic 
points to form a new set of points corresponding to panel and 
subpanel corner points. (Note, this is just a temporary set 
as the user may opt to use the A.P.R,. in the spanwise direction 
as well, in which case the section points do not necessarily line 
up with panel edges.) Subpanel points are always obtained by 
interpolation whether or not the A.P.R. has been activated. 

No matter which paneling option has been selected, basic 
points should be reasonably dense in regions of high curvature, 
such as near the wing leading ed'.,e. 

Several options have been provided for defining the basic 
points and these, jn combination with th( two-stage transforma- 
tion described .above, provide great flexibility when preparing 
the input. The options may be exercised at the section level 
so the input form may be changed from section to section. The 
options available at this time are described below and are con- 
trolled by the value of INPUT. INPUT values of 1 through 4 are 
illustrated in Figure 24 together with instances for their use. 

INPUT values of 1, 2 and 3 are used when a section lies in 
one of the reference planes of the chosen S.C.S.; in these 
cases we have a constant coordinate, x, y, or z, respectively. 
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Figure 24. Basic Point Input Options 1 Through 4. 





with one coordinate fixed we need input only two coordinates for 
each basic point, e.g., y and z when INPUT « i. provision is 
made to specify a third quantity to give a local adjustment to 
the "constant" coordinate, e.g., when using INPUT = 2 we may 
specify x, z and 6y. Usually the S-quantity is left blank 
(i.e., 0). The basic value for the constant coordinate is zero 
until the section points are transformed into the C.C.S., so 
the value of that coordinate in the C.C.S. must be provided in 
the transformation information. 

INPUT value of 4, v/hich requires all three components of 
each basic point position vector, is used when defining a com- 
pletely arbitrary section shape. 

INPUT values greater than 4 are provided and access a NACA 
four-digit equation to automatically generate a symmetrical 
section, the thickness/chord ratj.o having value INPUT/100 (i.e., 
there is a lower thickness limit of ,05c) . This option was 
originally used to chock out the geometry routines but it has 
been left installed as a possible future convenience--other 
equations could bo substituted easily. The coordinates are 
generated in the INPUT = 2 format, i.e., x, z with y = 0. 

Zero or negative INPUT values allow the present section's 
basic points to be copied over completely from any previously 
defined section. The section nunber is (-INPUT) except when 
INPUT = 0; the latter copies over the points from the section 
just completed. The section number specified is the absolute 
number from the beginning of the input and includes other copied 
sections as well as sections which may have been generated auto- 
matically. If the section counting becomes complicated, alter- 
native ways of copying are available as described later in 5.3.1 
The basic points are copied from the S.C.S. set (i.e., as origin 
ally specified) and are then transformed to the present C.C.S. 
according to the new section's transformation information. 
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5.2.3 Chordwise Regions 

The basic points defining a section may be assembled 
in a number of CHORDWISE REGIONS for the purpose of controlling 
the panel density and distribution on that section. In addition, 
the option on manual or automatic paneling is selected at the 
ehordwise region level, allowing the user to switch from one to 
another within each section vnerever he chooses. Chordwise 
regions are used only as an input convenience and are discarded 
in the program as soon as the surface paneling is complete. 

A chordwise region must end on a basic point called a NODE 
POINT, Figure 25. A NODE CARD, containing the chordwise region 
paneling information (see below) , inserted after a basic point 
in the input dock identifies that point as the end of a chord- 
wise region Node points are usually placed at "problem" areas 
where largr velocity gradients are expected to occur, e.g., flap 
hinge line,, loading edge, close-interference regions, but the 
user can place them wherever ho wishes to change from one panel 
scheme to another. Four typos of node point are provided at this 
time and are described below. 

The information on a NODE CARD consists of just three inte- 
gers. 

(i) NODEG identifies the node point and its type. 

(ii) NPANC is the numbc’ of panels to be generated by 

the A.P.R. in the chordwise region just com- 
pleted — a zero value gives manual paneling. 

(iii) ISPAC controls the form of rhe distributioii in the 
automatic paneling mode and is inactive 'in 
the manual paneling mode. 

(The C on the end of each quantity distinguishes the chord- 
wise from the corresponding spanwise quantities, which end in S, 
5.2.4. ) 
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Figure 25. Chordwise Regions on a Section. 


NODEC values of 1 or ?. specify the end of a chordwise 
region v/ith, resj, •’tively, continuous or discontinuous surface 
slope onto the nexu chordwise region. These values are, there- 
fore, used only on regions ending in the interior of a section. 
The last point on a section is specified by NODEC = 3 and is 
the only node point that must always be specified even if manual 
paneling has been selected. Negative NODEC values are also per- 
mitted and initiate a special copyi.ng routine described in 5.3.1. 

Four panel spacing options are provided in the A.P.R. The 
action of ISPAC values of 0, 1 and 2 is illustrated in Figure 26 
and is based on the cosine distribution giving increased panel 
density towards, respectively, the beginning and end, the begin- 
ning onjy, or the end only, of the region. Equal spacing 
throughout the re«jion is provided by ISPAC = 3. Coupled with 
the flexibility offered by the choice of chordwise region loca- 
tion. these spacing options have proved adequate so far; however, 
other options could easily be added should the need arise later, 
e.g., one based on increments in integrated surface curvature, or 
on increments in doublet value from a preliminary two-dimensional 
solution for the section. 

Clearly, node cards provide the user with an extremely 
versatile paneling tool. With one card deck of basic points 
defining the configuration geometry, he can, from run to run, 
change the form of the paneling simply by changing two integer 
values on each node card. Not only that, he can also move node 
cards within the deck (but not the node cards at section ends) 
or remove some or add new ones from run to run. This allows the 
user to concentrate his paneling in areas of interest, leaving 
other areas more sparsely paneled. It thereby provides a very 
effective use of the limited nuiober of panels available, yet, 
on a subsequent run a few small changes to the node cards allow 
the emphasis to be switched to another area without having to 
punch a new basic geometry c^rd deck. 
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s, - so -cos(e,))/z 

WHERE 8^ « (i - l)w/N AND N IS THE NUMIER OF INTERVALS REQUIRSO 

» ISPAC 0 



s, = SO - COS(0|)) 
WHERE 0j = (I - 1 )it/2N 

(b) ISPAC = 1 



Sj = S SSN(8|) 

WHERE 0| « (I ■■ l)n/2N 

(c) ISPAC = 2 


Figure 26. Spacing Options 0, 1 and 2 in the A.P.R. 


There is just one important ground rule for the use of node 
cards; the total number of panels (automatic and/or manual) on 
each section of a patch must be the same . The total is, in fact, 
the number of panel rows, NROW, for that patch. The program 
monitors the number of panels on each section and the calcula-^ 
tions are terminated with an error message should the user make 
a mistake. Provided this ground rule is satisfied, it is not 
necessary for the panel distribution to be the same from section 
to section — in other words, the number of chordwise regions and 
their node information can vary from section to section. The 
significance of this will be illustrated in 5.4.1. 

5.2.4 Spanwise Regions 

Sections defined within each patch may be assembled in 
a number of SPANWISE REGIONS for the purpose of controlling panel 
density and spacing in the spanwise direction. In forming span- 
vise regions, sections defined by the user take on a similar 
role to that of basic points in the chordwise regions. Although 
the options available for the spanwise regions are essentially 
the same as described for the chordwise regions in 5.2.3, the 
two are applied completely independently; for example, the user 
may request automatic paneling in the chordwise direction and 
manual in the spanwise direction. As in the case of chordwise 
regions, spanwise regions are u.sed only as an input convenience 
and are discarded once the paneling is complete. 

Spanwise regions must end at user-defined sections, called 
NODE SECTIONS, Figure 27. These usually coincide with kinks in 
tile spanwise direction on the patch planform, but the user can 
place one whenever he wishes to change the form of the paneling 
or to change between manual and automatic paneling in the span- 
wise direction. For convenience, the spanwise nr.ae information 
is included on the section card together with the section trans- 
formation information (5.2.2). The function of tht spanwise 
region node quantities, NODES, NPANS , ISPAS — distinguished from 
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INTERMEDIATE INPUT SECTION 



Figure 27. Spanwise Regions on a Patch. 


i 


63 


the corresponding chordwise quantities by ending in S — follov/s 
closely the description in 5.2.3, NODES, however, must be set 
to zero (blank) on the first section of a patch and on all inter- 
mediate input sections that are not node sections. (NPANS and 
ISPAS are then inactive.) The last section on a patch is identi- 
fied by a NODES value of 3, 4 or 5; 4 is used if the patch is 
the last one on a component and 5 is used if the patch is the 
last one on the configuration, in which case the present section 
completes the basic description of surface geometry. 

The total number of panels defined (manually or automatic- 
ally) across each patch in the spanwise direction is monitored 
by the program and becomes the number of panel columns, NCOL, for 
that patch. In view of the ease of generating panels, the code 
also monitors the running total of panels, and if a limit is 
exceeded, the calculation terminates with an appropriate error 
message. The limit is set internally by the storage capacity, 
but the user is given the opportunity to override that value 
with his own estimate of the total he intends to use for that 
case. In the event of an input error, this will avoid the in- 
advertent and expensive use of, say, 1,000 panels when the user 
intended using only 100. 


5 . 3 Special Routines 

The geometry routines described above may be applied for 
the complete configuration; however, special routines have been 
provided to reduce user input and, in particular, to avoid du- 
plicating information already supplied. These routines, which 
are described below, are optional. 

5.3.1 Copying Routine 

We have already seen (5.2.2) a copying facility acces- 
sible at the section input level. This copies over a complete 
section, including the chordwise region information, and has. 
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therefore, a rather limited application. More general copying 
routines are provided and are activated at the basic point level 
to copy STRINGS OP BASIC POINTS, rather than complete sections. 
This capability allows a new section to be assembled from parts 
of previously defined sections. Several strings of basic points 
may be assembled from a number of previously defined sections and 
the points selected need not follow the same direction as origin- 
ally specified. Furthermore, the copied strings of points may 
be intermixed with strings of manually input basic points to 
complete the new section. 

For this copying mode, the value of INPUT on the section 
card (5.2.2) must be in the range 1 to 4. The copying is acti- 
vated by inserting a NODE CARD having a NEGATIVE sign on NODEC. 
This is regarded as a DUMMY node card because it does not neces- 
sarily terminate a chordwise region (see below) . The negative 
value for NODEC determines the action at the end of the copied 
string of basic points. If NODEC = -1 or -2, then the last 
copied point becomes the end of a chordwise region on the new 
section and signifies, respectively, continuous or discontinuous 
slope onto the next chordwise region. We then continue to spec- 
ify further basic points, or, by inputting another negative node 
card, we can copy another string of basic points, and so on. If 
NODEC = -3, then the last copied point in the string completes 
the new section. 

If the user does not require a chordwise region to end 
at the last point in a copied string, then he sets NODEC = -4 
when he initiates the copy. When the string has been copied 
over, the program then expects to receive further basic points 
to complete the chordwise region or another negative node card 
can be used to copy another string of points, and so on. Clearly, 
if NODEC ~ -4, then the NPANC and ISPAC values on the NODE CARD 
are inactive and may be left blank. 

Whenever a negative node card is inserted, it must always 
be followed by a COPY CARD containing the following information 
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(four integers) defining the location of the required string of 
points, IPCH, ISEC, IB, LB. 

IPCH is the patch number containing the required points. 

ISEC is the section number relative to the start of that 

patch. 

IB, LB are, respectively, the first and last basic point 

numbers (inclusive) defining the string. The number- 
ing is relative to the start of the section ISEC. 

Thus, even in a complicated configuration, it is relatively 
easy to specify a string of basic points. 

This option offers not only an alternative to the earlier 
copying routine, but also a more general capability because the 
copying is initiated at the basic point input level, rather than 
at the section input le^^el. For example, the complete copied 
section need form only a part of the new section, it being pos- 
sible to have other basic points, both before and after the 
copied string. In addition to this, the ability to break the 
copying into strings of points allows a new distribution of 
chordwise regions to be selected. 

One restriction must be considered when using this copy 
routine — the new section's value for INPUT must coincide with 
the INPUT values on sections from which strings of points are 
to be copied. This restriction has not posed a problem so far, 
but if it does, it would not be too difficult to remove it. 

An example of the use of this copy routine will be shown in 
the test case in 7.1. 

5.3.2 Automatic Patch Generator 

Patches covering tip edges, flap edges, cutouts, etc. 
can be input by the user as ordinary patches, but this can get 
tedious. Optional automatic procedures have been installed which 
simplify this input by generating a complete patch within the 
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code. This AUTOMATIC PATCH GENERATOR, or A.P.G., is initiated 
at the patch input level by inserting a non- zero value for para- 
meter, MAKE, on the patch data card. The value of MAKE identi- 
fies the patch number on the edge of which a closing patch is to 
be generated. The sign of MAKE determines whether the new 
patch is on side 3 (positive) or side 1 (negative) of the basic 
patch. 

Consider, for example, a tip-edge patch. Here we have al- 
ready defined the patch representing the main surface. The end 
section of that patch provides the BASE SECTION from which the 
A.P.G. creates the new patch. Figure 28(a), according to user 
instructions, when the A.P.G. has been activated, the next 
card must contain the following: 

NPANC, ISPAC, KURV, NODES, NPANS, ISPAS. 

Referring to 5.2.3 and 5.2.4, the generated patch has one chord- 
wise region with NPANC panels spaced according to the value of 
ISPAC. It has one spanwise region with NPANS panels spaced 
according to ISPAS. The value of NODES must be either 3, 4 or 5, 
depending on the location of the patch in the input. The func- 
tion of KURV is described below. 

Sections defining the new patch are created automatically 
from the base section coordinates. The contour of each section 
generated may be either a straight line ("square-cut" tip) or an 
ellipse, depending on the value of the quantity, KURV, supplied 
by the user. If KURV is 0, sets of basic points are generated 
on straight lines joining upper and lower points on the base 
section. The same number of points is created even if the inter- 
val across the base section is zero (e.g., at the leading and 
trailing edges). Figure 28(a). 

If KURV is 1, the basic points are created on semicircles 
having diameter equal to the local "thickness" of the base sec- 
tion. 
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IF ICU*V>l(l.e. , X, , Y| ; 1=1, KURV) 


(b) PATCH WITH SEMICIPCULAR (KUtV - I) 

OR SEMI - ELLIPTICAL (KURV>1) SECTIONS 


Figure 28. Automatic Patch Generator. 


68 


If KURV is gi eater than 1, the basic points are created on 
semi-ellipses; the base section local thickness provides one 
axis, while the semi -axis is derived from additional user input. 
A planform shape is input using a set of coordinates, x^, 

1 = 1, KURV, defined in a convenient local coordinate systeift 
with origin at the trailing edge. Figure 28(b). The scale and 
point distribution are completely arbitrary, so the points may 
be conveniently measured from a planform view of the wing. The 
program scales the shape to fit the length of the basic section 
and interpolates to find the local semi-axis for each ellipse. 

KURV may also take a negative value. The generated patch 
then uses the copying routine (5.3.1) to locate the four sides 
of the new patch. In this case we input for each side a copy 
card containing the four integers, IPCH, ISEC, IB, LB, to locate 
four strings of basic points on previously defined sections. 

The A.P.G. then joins points on side 4 to points on side 2 using 
straight lines, i.e., sections on the new patch. The first and 
last sections on the new patch are taken directly from the 
strings of basic points for sides 1 and 3, respectively. I^ the 
copy cards for either of these sides has been left blank, then 
the appropriate first or last straight line joining sides 4 and 

2 become, respectively, the first or last section for the new 
patch. This option in the A.P.G. is useful for fitting in side 
openings left by flap-edge cutouis, Figure 28(c). Clearly, 
when KURV is negative, the value of the parameter, MAKE, is no 
longer important (as long as it is non- zero) ; however, the sign 
of MAKE still determines which way the patch is facing. Accord- 
ingly, we set MAKE = +1 or -1. 

One final point — the A.P.G. works directly in the component 
coordinate system (C.C.S.) even when copying strings of points. 
It always generates patches with positive IDENT (5.2.1). 
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5 V 4 Panels and S ubpanels 

5.4.1 Automatic Paneling Routine 

When the basic geometry has been specified, the panel 
ai'd subpanel corner points are assembled, patch by patch, A 
temporary set of chordwise points corresponding to panel and sub- 
panel corners is first assembled on each of the defined sections. 
This is performed in each chordwise region in turn (5.2.3) and 
interpolation is used when the A.P.R. is requested (i.e., when 
NPANC > 0) . Subpanel points are always generated by interpola- 
tion; if the A.P.R. has been requested, then the subpanel points 
are included as a set with the panel points. If manual paneling 
is being used, the subpanel spacing is based on normalized point 
subscript (Appendix A) ; this creates subpanel intervals more 
closely related to the changes in spacing in the user-specified 
panel points. 

The form of the interpolation used by the A.P.R. depends on 
the mimber of basic points available in the chordwise region, 
including the two end points. The code augments this number by 
taking a basic point from a neighboring chordwise region if con- 
tinuous slope has been specified onto that region (i.e., NODEC 
= 1 on i,his or the previous region) . The A.P.R. takes the 
available set of basic points and first eliminates zero length 
intervals, then, depending on thu mamber of basic points left, 
i.e., one, two, three or more, it uses, respectively, constant, 
linear, quadratic or biquadratic interpolation to generate the 
panel and subpanel points. 

When the temporary set of chordwise points is complete for 
all sections on a patch, corresponding points on each section 
are joined by lines called SPANWISE GENERATORS. The panel and 
subpanel points along each spanwise generator are then assembled 
in a similar way to that described for the chordwise direction, 
but now based on the spanwise region information. The interpola- 
tion routine is now applied along each spanwise generator in each 


spanwise region where the A.P.R. has been selected. The new set 
of (spanwxse) points are actual panel and subpanel corner points 
from which the panel and subpaiiel geometry is generated. 

The fact that we input just one set of spanwise region in- 
formation for a patch means that the same spanwise interpolation 
format is used on all the spanwise generators on that patch. 

Thus the A.P.R. in the spanwise direction has lost generality 
compared with the chordwise capability; however, this loss is not 
serious (and to avoid it would require considerably more input) . 

At this time, therefore, the combination of the chordwise and 
spanwise A.P.R. has the capability illustrated in Figure 29. 

The general character when using simple input sections (planar) 
is shown in Figure 29 (a) , while the availability of more general 
sections would allow spanwise stretching or compression in the 
paneling. Figure 29(b). 

5.4.2 Panel and Subpanel Geometry 

The four corner points, R^, i ~ 1, 4, specifying either 
a panel or subpanel quadrilateral are in the same sequence as the 
corners on the parent patch. Figure 30. From these points we con-, 
struct the two diagonal vectors 

Si = ^3 - ,4) 

R 2 = ^4 ■ ^'2 


The vector product of these diagonals produces a vector normal 
to the mean plane of the quadrilateral. We thereby construct 


the unit normal vector; 

n ' 


Si . Si/ls 


2 ' ;~.i - S 2 1 I dent 


( 6 ) 


The value of iDENT (5.2.1) is taken from the parent patch 
and ensures that the unit normal is always directed outwards 
from the surface into the flow field. The modulus of the 
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Figure 29. General Character of Paneling Offered by the A.P.R. 
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Figure 30. Panel or Subpanel Geometry. 
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diagonal vector product also provides the area of the quadrila- 
teral projected onto the mean plane; 

AREA = ^ 


The center point is taken as the mean of the four corner 
points ; 


Rc 




Ri ) /4 
" 1 


( 8 ) 


Two unit tangent vectors, m, are constructed, which, together 
with n form a right-handed orthogonal unit vector system for local 
coordinates. Thio system takes the center point, Rc, as origin. 

Tangent vector, m, is always directed from ^ to the mid- 
point of side 3 of the quadrilateral. Figure 30. (This always 
places m in the mean plane of the quadrilateral, even if the 
corner points are not co-planar.) Thus, 


m 




+ R^ /2 - Rc 




(9) 


With m and n known, we can construct 

^ = m ^ n (10) 


Next, the projections of the four corner points onto the 
mean plane are expressed in terms of the local coordinate system; 
that is, the relative position vector of from the local origin 
is 


E, 


= R, 


Rc 


( 11 ) 
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This has components, Ex^, Ey^, Ez^, say, in the local coordinate 
system, or E^^ == ^^lE' where Ex^ = E^ • I, etc. 

The quantity, Ez^^, is the projection distance of the corner 
point from the mean plane and indicates the amount of skew of 
the quadrilateral from its moan plane „ The magnitude of Ez^, 
which is the same for all four cornci points, should be kept 
small in relation to the size of the quadri].ateral . 

The quantities, Ex^, Ey^, i = 1, 4, define the flat projec- 
ted panel or subpanel in the local coordinate system and are 
used in the influence coefficient routine. 

Finally wo evaluate the half median lengths, SMp, SMq, for 
the quadrilateral. '^heso are (Figure 30): 


SMp ^ 
SMq - 


( ^2 f ^2) 

I £3 + E ^)/2 - EC 


( 12 ) 

(13) 


These are the half-lengths of the diagonals of the parallel- 
ogram which is always formed when the midpoints of adjacent sides 
of a quadrilateral are joined--even if the latter's corner points 
are not coplanar; the parallelogram lies in the mean plane of the 
quadrilateral and its area is half that of the projected quadri- 
lateral . 

Within each patch, the regular arrangement of panels and 
subpanels causes the adjacent side midpoints of neighboring sub- 
panels to coincide tixactly. This allows the SMp and SMq lengths 
of subpanols to be linked, respectively, in the chordwise and 
spanwise direction over the patch and thereby provides a c].ose 
approximation to surface distances between subpanel centers. 

The geometric data evaluated above (except for and D 2 ) 
are stored for each panel and its subpanels as a complete set. 

The arrangement of subpanols within each panel is always the 
same, see Figure 31. 
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5.4.3 Panel Neighbor Routine 

In order that a reasonable two-way interpolation and 
differentiation of the surface doublet distribution can be per- 
formed, it is important that we can quickly locate for each 
panel the set of four neighboring panels and their orientation. 
Each panel, therefore, keeps an array of four neighboring panels, 
NABORj^, i = 1,4 (in the same sequence as its sides), together 
with the adjacent side numbers of those panels, NABSID^, i = 1, 

4, Figure 32. Tne side number takes a negative value if the 
order of the side.'j on the- neighboring panel is reversed relative 
to the present panel — it is useful to regard this as a change in 
panel POLARITY. This reversal can occur when a neighboring panel 
is from a patch with a different IDENT (5.2.1), or when a panel 
takes a reflection of itself, e.g., at the plane of symmetry. 

Clearly, within the rectangular grid system of a patch lo- 
cating neighbors is easy; even so, the neighbor information is 
still stored to form a consistent system and to avoid repetitive 
calculation. Across the joints between patches, however, panel 
neighbors are not immediately available; for example, one panel 
may be neighbor to several smaller panels on an adjoining patch 
as we saw earlier in Figure 18. An automatic procedure has, 
therefore, been installed in the code which scans patch side 
panels in a search for possible i eighbors across patch joints. 

In this search only patches withxn the same assembly of compon- 
ents are considered. "Undesirable" neighbors are quickly elim- 
inated on the basis of relative geometry during the assembly of 
a short list of possible neignbors for each side panel. From 
this list of candidates one PREI ERRED NEIGHBOR is selected. 

At this time, preference is given to the panel whose control 
point lies closest to a normal plane constructed on the side 
panel. Figure 33. This plane contains the side panel's control 
point, unit normal vector and the side midpoint of the middle 
subpanel at the patch edge. 
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5.4.4 Subpanel Usage 

Two extremes of aubpanel usage are illustrated in 
Figure 32; on the left of the figure, panel size is fixed as 
subpanel density (i.e., niiinber of subpanels per panel) is in- 
creased, while on the right, subpanel size remains constant as 
the density increases. For a given surface area, the former 
case keeps the number of panels (and, hence, unknowns) constant 
as the subpanel representation becomes increasingly detailed, 
while on the right, the number of panels decreases rapidly with 
increasing subpanel density. In practice, we should fall be- 
tween the two, but the potential savings in the number of un- 
knowns, examined below, indicates we should perhaps lean towards 
the system on the right in Figure 34. 

For the purpose of evaluating the potential savings in the 
number of unknowns for the system on the right, we consider a 
surface represented by a 105 x 105 system of subpanels. We 
then assemble these into panels having 1x1, 3x3, 5x5 and 
7x7 subpanel arrays. Table 1 compares the corresponding 
number of unknovms with the number of influence coefficients 
(panel plus subpanel) evaluated for one velocity calculation. 

The latter assumes ten panels have their subpanel systems ac- 
cessed. The table also compares the approximate storage require- 
ment for the subpanel doublet multipliers (see 6.1) with that 
for the matrix of normal influence coefficients. The table 
shows the rapid reduction in the number of unknowns as the sub- 
panel density increases with the major reduction occurring for 
the 3x3 scheme, i.e., 89% compared with 96% for the 5x5 sys- 
tem. The corresponding reduction in number of influence coeffi- 
cient evaluations when going from the basic panels (i.e., the 
1x1 system) follovi’s a similar* trend? the 3x3 system offers 
an 88% reduction compared v/ith a 94, „ reduction for the 5x5 
.scheme. Going to higher subpanel densities than this decreases 
the benefits because of the high number of subpanels involved 
with the assumed 10-panel near-field set. 


81 


□ 


SmTANH DCWil fV 


• PANH Size 
COKtSTANT 

• NIiMKROr 
UNKNOWNS 
CONSTANT 

•DETAIL 

REPRESEN- 

TATION 

INCREASES 


*NUMM^ 

*SIS&KiATIOH 

CONSTAHT 




























(a) PANEL SIZE CONSTANT (b) SUIPANEL SIZE CONSTANT 


Figure 34. Extremes of Subpanel Usage. 
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The additional storage neede.d for the subpanel doublet mul- 
tiplier is clearly insignificant in comparison with the savings 
in storage for the matrix of influence coefficients. 

These comparisons are based on a simplified arrangement and 
should be regarded only as a rough guide. Even so, there is an 
obvious attraction to use a 3 x 3 system and little point in 
going to 7 x 7 densities or higher. For this reason the new 
coding has an upper limit at the 5x5 subpanel scheme. 

5. 5 Wake Routines 

5.5.1 Initial Wake Geometry 

Wakes are formed after all the surface patches have 
been paneled and the neighbor information stored. The user 
identifies strings of WAKE-SHEDDING PANELS, the side geometry 
of which defines the FIXED BASELINE of each wake. At the end of 
each string of wake-shedding panels, the user has the option 
of defining the initial (i.e., prior to wake relaxation) stream- 
wise geometry of a line on the near-wake (4.2) using a set of 
BASIC WAKE POINTS. The function of these points is similar to 
that of basic points defining chordwise sections on patches 
(5.2.2) . Node cards are used here also and allow the user to 
select wake panel density and form of distribution in accordance 
with the expected location of the relaxed wake; in this way, 
wake panel detail mny be used efficiently in relation to the 
expected wake curvature and surface interference . As a minimum, 
one basic wake point must be provided on each streamwise line 
defined — this point corresponds to the downstream end of the 
near-wake. (The upstream end is taken automatically from the 
fixed base line.) Multiple basic wake points are essential only 
in the case of multiple high-lift devices and allow a represen- 
tative initial wake geometry to be defined which should reduce 
the number of wake shape iterations later in the calculations . 
Figure 35 illustrates the case of a v/ing with slat and slotted 
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igure 35. Basic Wake 



flap and shows a typical set of basic wake points. 

When the basic wake information has been supplied, the pro- 
gram generates streamwise sets of wake panel corner points ac- 
cording to the user's instructions on the node cards. A set of 
points is generated corresponding to each wake-shedding panel 
corner point on the fii’od baseline and ends at the downstream 
edge of the near wake. Linear interpolation is used in both 
streamwise (between basic points) and spanwise (between stream- 
wise lines) directions for the purpose of generating these 
initial sots of wake panel corner points. The user should boar 
this in mind when selecting the number and location of both 
basic points and streamwise lines. It must be emphasiKod, how- 
ever, that this information is used only to define the prolimin' 
ary wake for the purpose of the first solut ioiv-- thereafter , the 
wake relaxation routine v;ill redefine the wake geometry at each 
iteration. 


5.3.2 W aicG Pa iw 1 s and S u t)p a n^e 1 s 

Thi‘ program processes the streamwise sets of wake panel 
corner i)oiuto and generates subpanol corner points using two- 
way biquadratic interpolation. V»'ake panel and subi^anel parameters 
are then formed as in the case of surface patches (5.4.2). Al- 
though there is an obviovis oimila ’ity between a patch and a wake, 
the doublet distribution on the lattor is loss complicated as it 
is constant along the (streamwise) columns. Doublet multipliers 
(see later in 6.1) associated with v;ake subpanels are therefore 
dependent only on the spanwise geometry, which change.o with each 
wait' relaxation iteration. 

The doublet value on each column of v/ake panels is the dif- 
ference between the values on the corresponding wake- shedding 
pajiol and its neighbor across the shedding line. This neighbor 
relationship across the shedding line is terminated once the 
wake has been formed, the doublet distribution then passes 
smoothly onto the wako fx'om both sides of the shedding lino. 
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5.5.3 Wake Relaxation 

After the initial processing of the input/ each near- 
wake shape is alwa* a stored in the form of the streeunwise sets 
of wake panel corner points. After a doublet solution has been 
obtained, velocities are calculated at each of these points, 
Figure 36. Each set of points is then relaxed into the local 
flow direction as defined by the local velocity vectors. These 
calculations proceed from the (fixed) first point in each set 
and finish at the downstream end of the near-wake. 

Based on the new sets of points, wake panels and subpanels 
are regenerated after each wake relaxation. The wake doublet 
values are assumed to move with the center subpanel in each wake 
panel. In this way the stretching and contraction of the wake 
affects the vorticity level when the spanwise gradient of the 
doublet distribution is evaluated. 


C “ 
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6.0 SURFACE DOUBLET DISTRIBUTION 

6.1 Doublet Multipllera 

In order to evaluate the doublet values at subpanel centers, 
we form a two-way biquadratic interpolation through neighboring 
panel doublet values. The interpolation is based on surface 
distances, p, q, respectively in the chordwise and spanwise 
directions over each patch. Surface distances are normalised 
to a, 3, respectively, using local intervals between control 
points. Figure 37(a). Intervals in p and q are evaluated on 
the basis of straight-line distance across subpanels , i.e., by 
using the stored SMq information (5.4.2). The two-way bi- 

quadratic scheme requires a 4 x 4 set of panels for interpola- 
tion within the central area, Figure 37(a). One 4x4 set, 
however, only covers one quarter of the panel; four such sots 
are, therefore, required to cover all the subpanels on a panel. 
Since the forr sets overlap they can be selected from a 5 x 5 
panel set having the panel at the center, Figure 37(b), i.e., 
within this 25-panel set, the panel itself is always at location 
13. with this arrangement, the selection of panel sets for sub- 
panels within each panel follows a common rule throughout the 
configuration. Not all the subpanels require a 4 x 4 panel 
set to define their doublet value; a subpanel which lies on a 
grid line joining panel control points requires only 4 panels; 
for example, with the 3x3 subpanel scheme shown in Figure 
37(b), subpanel number 3 is related to panels at locations '3, 

8, 13, and 18. In particular, the middle subpanel takes the 
panel doublet value so it has just one multiplier which has unit 
value. 

The doublet value at location a, 3 is given by 


y (a, 3) 


* E 

i»l,4 
j*lf 4 


y 


ij 



(14) 
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Figure 37, 


Two-Way Biquadratic Interpolation for Subpanel 
Doublet Values. 
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where i and j refer, respectively, to the row and column numbers 
in the 4x4 panel set and are the doublet multipliers evalu- 
ated at a, p and are given by 

°ij ■ '^ij (15) 


(This is a simplified scheme which may give distortions in ex- 
tremely irregular grid patterns— an alternative, more rigorous 
(but also more complicated) scheme is held in reserve should the 
simple scheme fail,) 

The d^j and e^^^ are the column-wise and row-wise biquadratic 
multipliers, respectively, evaluated as follows: 


first define 

“1 j " ^Pl Ap 2 . 

3 / 3 

(16) 


“ 2 j “ ^^^ 2 . / "^^ 2 . 

(17) 


where Ap^j = are the surface length increments; 

hence, 

d^. « Gl(a,a^ ); d 2 * = G2(a,a, ,«2 ) »' (18) 

j j j 

*^31 “ <^2(1 - a, 1 - tt 2 » 1 - a, ); (19) 

j j 


^ 4 j = Gl(l - a, 1 " 0 I 2 ) ; (20) 

with G1 and G2 being the general biquadratic multipliers (Appen- 
dix A) . Similarly, in the 3 direction 


( 21 ) 

5i2 = (Alj 2 + * 1 ^ 3 ) j ( 22 ) 
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Hence , 


e^3 - G2(l - a, 1 - a^2' ^ 

®i4 * ^ " ^i2^ • (25) 


At the beginning of the calculation, the panel doublet 
values, (solution part) are unknown and so the subpanel 
doublet multipliers, j > must be stored to be later applied 
to subpanel influence coefficients (6.3.3). 

6.2 Augmented Patch 

To facilitate the collection of the 25-panel sets in a 
regular manner without problems at patch edges, an AUGMENTED 
PATCH is temporarily formed for each patch in turn. The aug- 
mented patch has a two-panel deep fringes of panels surrounding 
the basic patch. Figure 38. The fringe panels are assembled 
from the neighboring panel information and the complete set of 
panels and grid distances, p, q, are formed for the entire 
augmented patch. Each panel in the basic patch can then quickly 
locate its 25-panel set together with the grid distances, even 
if it is a one-panel patch. 

In forming grid line distances within the fringe areas, 
the orientation of the neighboring panels must be considered. 
Four panel sides and the panel polarity (4.5.3) lead to eight 
possible orientations of a neighboring panel across a patch 
side. All these possibilities are covered by the coding, but, 
further refinement is needed when evaluating surface distances 
in the fringe area in cases where there is a large mismatch in 
the neighboring panel alignments. 


92 


PANH KQUCHCE IN 



(DO NOT NECESSAMLY LINE 
Uf WITH lASIC PATCH GRID) 


Figure 38, Augmented Patch Scheme for Doublet Interpolation. 
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6.3 Influence Coefficient 

The influence coefficient routine for the doublet singular- 
ity model la arranged in three parts for application in far-field, 
middle-field or near-field velocity calculations. So far, the 
distances to the boundary of each of these regions has not been 
fully explored, but approximate guidelines are indicated below. 
These distances, which are measured from the panel control point, 
are given in terms of panel "size"; this is defined as the square 
root of the panel area at this time. 

6.3.1 Far-Field 

In the far-field, say beyond ten panel sizes away, 
the doublet distribution on a panel is regarded as piecewise 
constant. The influence coefficient is then that of a quadrila- 
teral vortex defined by the panel's four corner points. The 
vortex strength, F, is the same as the panel doublet value. 

The velocity induced by each of the four sides of the quadri- 
lateral vortex is obtained from the Biot Savart law; 

a ^ b(a + b)/{&b(ab + a • b) } (26) 

, etc. 

where a and b are the position vectors of the velocity calcula- 
tion point relative to the start, and, respectively, of the 
straight line segment. The lengths of these vectors are deno- 
tated by a and b, respectively. 

The form of the induced velocity expression, which was 
developed during the course of this work (Appendix F) , eliminates 
the numerical problem associated with earlier forms when calcu- 
lating a velocity close to the extended line of the vortex. The 
new expression passes correctly through zero without special 
treatment for this condition. 


V = 


T_ 

4tt 


where 


a = 


a 
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6.3,2 Middle Field 

Between r say, three and ten panel sizes away, a panel 
is regarded as having a two-way linear vorticity distribution. 
The vorticity value and slope are derived from the doublet dis- 
tribution by passing quadratic curves in the p and q directions 
through the neighboring panel doublet values. Figure 39(a). 

Using the nomenclature shovm in Figure 39(a), the gradient 
and second derivative of the quadratic curves in the p and q 
directions are evaluated at the p 2 mel center as follows: 


- Md ^Pb/' 

(aPoCAPo + APb>} 

(27) 

^ 2 { <>^d/'^Pd 

Pb/''Pb> / <*Pd *Pb' 


” Wfl / 

APb)| 

(28) 

aq = {“C - U„(Aq;^ - Ag^) ^ 

/ Ma } / Aqc 

Aqc / ] 

(Ma'a^a Me)} 

(29) 


^ - 2 I /iiq^ + p,, j + 4q^,) 

- I (Aq^ Aq^) } 
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Figure 39. Evaluation of Vorticity Value and Gradient on a Panel. 
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Th« vortlcity coniponentB In the panel's local coordinate 
system (see Figure 39(b)) are obtained at the panel center: 


and 


Yx “ Iq 


Y'x 


IhL 

3q^ 


(Ji) 


(y gradient in the q direction) 

(32) 


■( m • t •“ 

\ 3q - - 3P 


}/i • t 


(33) 


and 


/iin ra . t - 


ra • t - \ I a • ^ 

3q=^ 3p^ '/ 


' tirr 2 


(34) 


Where t is he unit vector in the direction from the panel 
center to the midpoint of side 2. 

The linear vorticity influence coefficient for the panel 
is accumulated by considering each side in turn, applying the 
model shown in Figure 40; for simplicity, this Illustration 
and the following description are for the case of y , but the 
Yy value is treated in a similar way. 

The basic linear vorticity model has two serai-infinite 
strips of opposing vorticity separated by a swept line (i.e., 
a panel side) . The strips are aligned with the vorticity vector 
while the vorticity gradient is normal to that direction (Figure 
40) . 

The velocity induced by the strip is 


V 


4it 


I vm + wn I 


(35) 
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where 


V 




T + y\ a^iL 


XJ) 


and 


w » Yxl^^ “ “ ^'x I + X (b " a)/e^ + a,^T) 


T » tan 


-1 


®n‘“Pb - ‘‘Pa'/'PaPb + "n' **=' | 


J « 


^yFl'^ - \ - *=)i>/V^ » - ®m - **n> 


L *= 


n 


(b - b^) (a + aj^)/(b + bj^)/(a 


a^) 


Pa = ^n'^ “ 


Pb = X(b^2 f a^2) - b^ b^ 


"'^xl ~ ^3c ^ x^^a ^ 


03 = 1 + X^ 


^1 “ " ^ 5 , 
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X » (Xj^ - 3c^)/yj3 “ Y^) slope of joint line 

a. , a , a , bo , b , b are the components of position 
fcm nA»mn 

vectors at a and b (Figure 40) in the panel's local coordinate 
system. 

Applying this model to all the sides of a (closed) polygon 
causes the vorticity to reinforce (double) in the interior re- 
gion and to cancel everywhere else. (The factor of 2 has been 
taken into account in the influence expression above.) 

6.3.3 Near-Field 

For velocity calculations within about three panel 
sizes from a panel's control point, the panel's subpanel set is 
accessed. (This includes the case for the influence of the 
panel on its own control point.) Each subpanel uses the linear 
vorticity model described in 6.3.2, except that here the vorti- 
city value and slope are evaluated using neighboring subpanel 
doublet values rather than the panel values indicated in Figure 
39. 

As each subpanel's influence coefficient is evaluated, it 
is fac ,ored by its set of doublet multipliers (6.1) to get the 
corresponding contributions for the associated local panel set. 
These contributions are then accumulated in the matrix of in- 
fluence coefficients. 


7.0 TEST CASES 

A number of preliminary tests have been carried out to 
check the working of the modified routines. wo of these cases 
are described below » but further tests need to be carried out, 
especially for general configurations. In particular, experi- 
mental data in part-span, high-lift devices is needed for ccm- 
parison purposes to thoroughly check the code. 

7.1 Geometry Code 

As a preliminary check of the paneling capability, a general 
configuration was assembled having part-span, high-lift devices. 
Figure 41(a). Calculations for this case were terminated after 
the panels had been generated. The configuration consists of 
three separate COMrONENmf? (limit is 10) represented by twelve 
PATCHES (limit is 100) . 

The slat is represented by three patches. Patch 1 covers 
the main slat surface — trailing edge through leading edge back 
to trailing edge, and root to tip. It is defined by two chord- 
wise sections — one at the root, and one at the tip. The tip- 
edge Patches, 2 and 3, are generated automatically (5 '*-2) from 
the root and tip edges, respectively, of Patch 1. Thet<i two 
patches were specified to be flat. The option provided to gen- 
erate a streamlined tin-edge patch is exercised on Patch 7 — 
the wing tip; here, an automatic tip-edge patch was requested on 
the outboard edge of Panel 6, but the planform contour was de- 
scribed. The contour description is a set of points (any spac- 
ing and any scale) going from the trailing edge to the leading 
edge. The program generates tip-edge sections using semi-el- 
lipses based on the local thickness of the basic patch edge 
(i.e.. Patch 6), and the local offset of the tip contour plan- 
form (after internal scaling and interpolation) . 

Other flat-edge patch options are exercised on Patches 8, 

9, 11, and 12. For patches 8 and 9, the partial section copying 
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Figure 41. Tests for the Geometry Code. 


option (5.3.2) was used; this special input requires 4 sets 
(one for each side of the now patch) of 4 integers. The four 
integers include the patch number, the section number on that 
patch, and the first and last point subscripts on that section 
which are to be copied for the new patch. A similar option is 
employed to define the first section of Patch 5: here, the 

majority of the section has already been defined on the last 
section of Patch 4. The f: st section of Patch 5 can, therefore, 
be defined by reading in : ' ivt points (seven in this case) to 

define the cove contour it»«e the detail in Figure 41) followed 
by a copy (5.3.1) of the set of common points from Patch 4. 

Figures 41(b) through 41(e) show plots of the panels and 
subpanels generated. Figure 41(b) shows the general view of 
the panels on the complete configuration and indicates the di- 
rections of detailed views shown in Figures 41(c), (d) and (e) . 
Figure 41(c) gives the detail of panels and subpanels on the 
wing tip. Patch 7. The closer geometric representation offered 
by the subpanel scheme is obvious in this case. (Note: the 

control point conditions for each panel are taken from the 
central subpanel on that panel. Also, the panel influence on 
itself always uses its basic subpanel set.) 

Figures 41(d) and 41(e) show similar details for Patches 
10 and 11 and for Patch 8, respectively. Patch 9, at the out- 
board end of the flap cutout, is very similar to Patch 8. 


7 . 2 Wing Case 

The modified potential flow code was applied to a rectangu- 
lar wing cf aspect ratio 2 to check the routines through to the 
pressure calculation. The wing section was the 11.1% t/c, 

Boeing Section TR 17 and angle of attack 5.73°. Figure 42(a) 
shows the chordwise pressure distribution at .125 semispan cal- 
culated using panels distributed. in a 24 x 4 array on the main 
surface patch and a 2 x 12 array on a tip patch with semicircular 
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Figure 41. Continued. 
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sections. For comparison^ Figure 42(a) also shows solutions 
from the original VIP3D program (36 x 4 array) and also from 
the USSAERO program (Ref. 9). There Is very close agreement 
between all three programs. This Is very encouraging because 
the doublet solution used a less dense panel system than the 
others . 

The tip patch paneling In the doublet model allows pres- 
sures to be calculated round the tip edge. Figure 42(b) shows 
pressure distributions plotted In the spanwlse direction from 
lower surface round the tip and back along the upper surface 
at x-wlse stations .0036 and .889. Values are plotted from two 
panel distributions, one with 4 equal spanwlse Intervals and 
one with 6 spanwlse Intervals with cosine distribution giving 
Increased density towards the tip. The latter Improves the 
matching In panel size between the main surface and tip patch 
compared with the first case which has panel size ratios of the 
order 50 passing onto the tip patch? this probably accounts for 
the discrepancies between the two soJutions near the tip in 
Figure 42(b). Large and sudden changes in interval size can 
cause numerical eri.or when inteirpolating or differentiating the 
surface doublet distribution. 

At the forward station, the spanwlse flow from lower sur- 
face onto upper surface clearly has a monotonically decreasing 
pressure. Towards the trailing edge, however, the upper surface 
suction level has disappeared while a peak suction has developed 
on the tip surface. Figure 42(b). At this station, therefore, 
the spanwlse flow is suddenly faced with a strong adverse pres- 
sure gradient as it climbs round the tip edge and will lead to 
the conditions for tip edge separation. 
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8.0 CONCLUSIONS AND RECOMMENDATIONS 

An investigation in two-dimensional flow demonstrated that 
a doublet subpanel technique has the required behavior for cal- 
culating interference pressures in a vortex/surface interaction 
caset the close-approach problems associated with earlier panel 
methods is essentially removed without increasing the panel den- 
sity. Because the subpanel model provides a better representa- 
tion of curved surfaces and smooth singularity distribution r it 
gives the effect of increased panel density without increasing 
the number of unknowns. Subpanels are accessed only on panels 
close to the velocity calculation point. 

Results from preliminary test cases of the three-dimensional 
form of the doublet subpanel technique have been encouraging. 

The VIP3D geometry routines have been extensively modified to be 
compatible with the technique and also to allow application to 
general high-lift configurations. For this purpose a .versatile, 
user-oriented scheme has been developed based on multiple patches 
of panels. Preliminary test calculations have shown very close 
agreement with solutions fror.' the original VIP3D singularity 
model. 

Further test cases need to be performed, but experimental 
measurements on general high- lift configurations are required 
having detailed pressure distributions and flow visualization. 

In particular, cases having part-span, high-lift devices need 
to be examined. 
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10.0 APPENDICES 


A. Biquadratic X.iterpol ation 

The biquadratic interpolation scheme described below is 
applied in a number of routines in the new code. Its simple 
multiplier form is very convenient to use and yet it is a "con- 
strained '* cubic, i.e., it cannot oscillate wildly. Experience 
with the routine over a number of years in the codes of Refer- 
ences 2 through 7 have shown it to be a relizdble method. 

Given a set of position vectors, Pn, n»l,2,...,N, defining 
a smooth space curve, we wish to interpolate for additional 
values in, say, the interval between P^ and Figure Al. 

We first generate the integrated contour length, 
point from the beginning of the curve, i.e., from P^^. For con- 
venience, the straight segment length is used across each inter- 
val, i.e., s =* I P . 1 " P I / but arc lengths could easily be 


s^, to each 


n 


substituted — or indeed, any other parameter that varies smoothly 
and monotonically along the curve, i.e., without introducing 
multiple value problems. For example, the point subscript in- 
terval is used in some parts of the program. 

Next, we generate two quadratic curves: q, (a, a, ) passing 

through points P » ?-n+l' 3,2 ^ P^issing thrcugh 

points Pjj 4 i» and ^^gure Al. 

The normalized i, erpolation parameter, a, ranges from 0 
to 1 in the n^l interval, and has value 


a 


1 at p , 
In £n-l 


and a« at P„,« 
2n -n+2 
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Figure Al. Biquadratic Interpolation. 


where 


a 


a 



(s 


I.-1 - V/'^n+l - 


n 


and 


*n+2 




Sn) . 


In the n^^ interval, we take a linear combination of 
and ^2 define the biquadratic interpolation curve there: 

r»3t) * 


n 


The biquadratic is, therefore, a cubic, but it is constrained to 
lie between two quadratic curves. It can't therefore, behave 
wildly. 

The value of a for a point distance, s, from the start of 
the curve (but located in the interval) is 


a 


s/(s 


n+1 


“ s ) 
n 


The form of the interpolation curve can be expressed in 
terms of biquadratic multipliers, Gl, G2, applied to the four 
local position vectors: 


P(a) == Gl(a,a^ ) + G2(a,ttj^ , o >2 ) 

n "^n n 


+ P^^j^ G2(l - a, 1 - U 2 f 1 ” ) 

n n 


+ P ^^2 Gl(l - a, 1 “ a, ) . 


n 
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The forms of 61, 62 are; 


61 (a, b) « a(l - a)V^b(l - b) } 

62(a,b,c) » (1 - a) [l - a(l - a)/b - a*/c} , 

These multipliers, based on the linear combination of tv;o 
quadratics, give continuous slope and a piecewise linear — but 
not necessarily continuous— variation of second derivative across 
each interval. A similar set of multipliers has been formed 
which gives continuous second derivatives, but it has not been 
thoroughly checked out at this time. 

The 6 multipliers can be differentiated to give the tan- 
gent vector; 

t(«) “ £n-l ) + £n “l ' “2 ^ 

■‘‘n ■‘"n ^n 

- H2(l - a, 1 - , 1 - ) 

n n 

— n +2 ^ °^2 ^ 


Where the tangent multipliers are: 

Hl(a,b) = (1 - 4a + 3aM/{b(l - b)} 

H2(a,b,c) = (4a - 3a^ - l)/b + (3a^ - 2a)/c - 1 
Note; this is not a unit vector; lt(a)l = 3 — . 
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The G multipliers can also be integrated to give the area 
under the curve between the n^^ and (n + 1)^^ points: 

in ■ [ En-1 > + £n ' “2 > 

L n „ n 

n 

+ P2(l - a , 1 - 0.1 ) + PIU - ">2 )] i 

n _ n 


where the integral multipliers are: 

Fl(b) « l/{12b(l - b)} 

E’2(b,c) « {1 - (1/b + l/c)/6}/2 

ds 

Thxs assumes that the value of is constant over the interval. 

In this case, therefore, it would be an advantage to use the arc 
length intervals rather than straight line intervals when cal- 
culating surface distances. 

The P multipliers are used in the Integration of surface 
pressures to obtain forces and moments. 
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B. Numerical Integration for Doublet Influence Coafficient 

A doublet singularity can be regarded as the limiting condi 
tion when a source and sink of equal strength became coincident 
in such a way that the product of their strength and distance 
apart remains constant. This constant is the doublet strength, 
y, and the unit vector, n, in the direction source-to-sink de- 
fines the axis of the doublet. The velocity induced by the point 
doublet is then 

( 3a • n a \ 

V « -Ji_ n - = 

4ira® ( a^ ) 

Where a is the position vector of the velocity calculation 
point relative to the doublet. 

The velocity induced by a doublet distribution on a flat 

#«« « 4S <n> M A Ct •• T1) *1 M X. 1 m, mm— 

ouxxaCcf £x^ui.c3 OJL f iiao lane 


1 3a(a,3) * k a(a,3)) u(a, 3 ) dpda 

Y - ■kjjh- — * 


(a, 3) 


a® (a, 3) 


where 


and 


a^ (a, 3) “ (a, 3) + z‘ 


r^(a,3) = (x - a)^ + (y -• 3) 


The component of velocity noarmal to the surface has polar 
symmetry about the doublet axis, k, and the integrand is 


Iw (r) 


r^ - 2z^ 

{j.2 + z2}5/2 
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This has the form illustrated in Figure B2 (a) . A number o£ 
interesting features can be evaluated. 

The negative peak value occurring at r « 0 depends on the 
height, z, of the point zd)ove the sheet. 

Iw(o) w ^ 
z* 

The radius at which Iw passes through zero is 
r^^ " /2|z| 


The secondary (positive) peak value is 


Iw(r2) 


55/2,. 


and is located at 


^2 - 2 z| 


It is important to note that both Iw(o) and IW(r 2 ) go singular 
(in opposite directions) as z tends to zero. 

The velocity component tangential to the sheet also has 
polar symmetry and is radial. The integrand for this has the 
form 


Iv(r) 




Which is illustrated in Figiire B2 (b) . 
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The peak value is 


IvCig) 


48 


and occurs at 


rg « |z|/2. 


The locations of the peaks in both velocity components 
were used to define separate regions in which the Romberg inte 
gration technique, Reference 10, was applied. Transformations 
were applied in each region according to the local behavior of 
each integrand. The aim was to obtain a more linear variation 
across each region. This would allow the Romberg Integration 
technique to converge more rapidly since it is based on the 
trapezoidal rule. 
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C. Linear Vorticity Influence Coefficient in Two 
biroensidna 

There are a number of two-dimenaional potential flow codes 
based on linear vorticity panels; however, the influence coeffi- 
cient used in the present work was formulated differently from 
earlier forms and is based on the panel center rather than the 
panel ends, Figure Cl. This is more convenient vjhen working 
with the surface doublet model. 

The induced velocity is 

V « t + Vjj n 


Where the velocity components tangential and normal to the 
panel are, respectively: 

''t ■ W Yo ^ 

''n “ ^ L + Y'(s - T + I. 

where Yq and y* are, respectively, the vorticity value and 
gradient at the panel center, and 

T = Tan”^ a^s/(a^ - ,25s^) 

L * I (a" - a^s + .25s^)/(a^ + a^s + .25s^) | 

a^ “ a • t 

a = a • n 


and s is the panel length. 
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D. Neuniann Boundary Condition Appliad to a Surface Doublet 

H55ir 

When the Neumann boundary condition ia applied to a aurface 
doublet diatribution, the aolution can diverge near the airfoil 
trailing edge. To find the aource of the problem, we exzunine 
the aimplified situation where the doublet sheets on the upper 
and lower surfaces are close to each other and are parallel, 
Figure Dl. Consider the Neumann boundary condition applied to 
a control point on the upper surface. The major terras are 

^(Uu ” * Vjj 


Where are the local upper and lower doublet values, 

respectively, and is the local normal component of the onset 
flow. 

Next, divide the doublet values into siTometrical and anti- 
symmetrical parts: 


and 





y 


A 


The equation above then becomes 


‘“{“s "a - - V { ■= \ 


or 


u = V 

N 
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Figure Dl. Parallel Doublet Sheets in Close Proximity. 
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J 


That is, the syiomatrical part, J4„, has disappeared. In 

B 

this simplified situation, therefore, the symmetrical part 
could take any value and yet the Neumann boundary condition would 
still be satisfied. Thus is indeterminate. 

In practical cases, this condition is rarely met with ex- 
actly; however, it is approached in the neighborhood of sharp, 
trailing edges. In these circiomstances the Neumann boundary 
condition is weak, and we sometimes observe that the upper and 
lower doublet values deviate from the exact solution, but they 
move together. The difference between them, i.e., the antisyra- 
metrical part, is generally close to the exact value. This 
numerical ill-conditioning is especially serious when using 
iterative solution techniques. 
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E. Streamline Calculation Routine 

A numerical procedure was developed for calculating 
streamline paths. The procedure is based on finite intervals, 
a mean velocity vector being calculated in the middle of each 
interval as we proceed along the path, Figure El. The velocity 
calculation point, RP, is obtained by extrapolation from the two 
previous intervals on the*Lnsis of constant rate of change in 
the velocity vector direction,* i.e., 

^ + .5 t 


where R^^ is the previous point calculated on the streamline, 
is the present interval length, and t is a projected tangent 
vector. 


t = 



V 

'-n-2 ) 


+ t 


n-1 


where t^_ 2 ^ = -Ty — — r t etc. Yn-i velocity calculated 

|— n-1 

in the middle of the previous interval. 

RP does not necessarily lie on the streamline path. 

When we have calculated the velocity, V , at RP we use this to 
evaluate the next point on the streamline, i.e.. 


R 


n+1 


= R 


-n 


S V 
n — n 

|2nl 


An automatic procedure is included in the routine which 
changes the interval length, S^, in accordance with the cnange in 
tangent direction. If a large change in direction is calculated, 
the value of is decreased and the calculation is repeated 
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GIVES LOCAL TANGENT VECTOR 


Figure El. Streamline Calculation. 


until the change in direction is within a specified amount. 

If the change in direction is smaller than a certain amount, then 
the interval length for the next step is increased. 
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F. Vortex Segment Influence Coefficient 

The fiimiliar expression for the velocity induded at a 
point, P, by a straight vortex segment is, referring to Figure 
PI; 


V 



cos 0^ 


+ cos 02 



n is the unit vector normal to the plane containing the segment 
and the point, P. 

A more convenient form for three-dimensional analyses 
(Ref. 4) is 


tr = 


JL 

4it 


(a 


b) 




/i 


ab / 


which avoids the evaluation of trigonometric quantities. 

Both expressions have a numerical problem when P approaches 
extension of the segment. (The case where P approaches the 
segment is a separate problem requiring special treatment, such 
as a core model . ) ^ 

The computer is then faceci with dividing one small number 
by another small number. The result should pass smoothly through 
zero as P passes through the in-line condition. In practice, 
round-off in the computer creates spurious results which neces- 
sitate a special local treatment. 

A solution to this problem has been found in the present 
work by rearrangement of the second expression. This can be 
written 

j, a ^ b(a + b) (1 - cos 0 ) 
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The cause of the numerical problem is clearly the sin^O 
term in the denominator — 0 being 0 when P is in line with the 
segment. This term can be cancelled after multiplying numerator 
and denominator by (1 + cos 0), leaving; 

a ^ b(a + b) 
a^b^' (1 + cos 0) 


a A b(a + b) 
ab (ab + a * h) 


or 


V 


r 

I? 


V = 


T_ 

4u 


This expression passes through zero correctly without special 
treatment as P passes through the in-line condition. 

The corresponding form for the semi-infinite vortex is 

V = - — ^ i 

~ a (a + a • t) 


where t is the unit vector along the vortex. 
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